Library Loading for the analysis

library(data.table)
library(tidyverse)
library(rtracklayer)
library(DESeq2)
library(Biostrings)
library(Rsubread)
library(pheatmap)
library(gridExtra)
library(ggrepel)
library(CLIPanalyze)
library(ggplot2)
library(RColorBrewer)
library(mixOmics)
library(plotly)

Loading merged data

Load the rda file from merging all peaks and filtered peaks from HF/HFK individual CLIPanalyze call. Counting was re-do using featureCount and DESeq2 were called to do differential analysis.

dir.create("PDF_figure/CLIP_DESeq_12232019_merged", showWarnings = FALSE)

load("../Merged_Analysis/merged_peak_analysis.rda")

plotMA(dds.genes,
       main = "MA plot for all HEAP vs IC\n in genes outside peaks")

pdf("PDF_figure/CLIP_DESeq_12232019_merged/MAPlot_all_HEAPvIC_OutsidePeaks.pdf",
    height = 4,
    width = 5)
plotMA(dds.genes,
       main = "MA plot for all HEAP vs IC\n in genes outside peaks")
dev.off()
## quartz_off_screen 
##                 2
plotMA(dds.peaks,
       main = "MA plot for all HEAP vs IC in peaks,\none-sided test")

pdf("PDF_figure/CLIP_DESeq_12232019_merged/MAPlot_all_HEAPvIC_InPeaks.pdf",
    height = 4,
    width = 5)
plotMA(dds.peaks,
       main = "MA plot for all HEAP vs IC in peaks,\none-sided test")
dev.off()
## quartz_off_screen 
##                 2

Peak filtering

No peak filtering was done here as peak was filtred and selected in their individual analysis.

HF/HFK comparison

plotMA(hfk.hf.res,
       main = "MA plot for HFK over HF\nin all peaks (significant over input)",
       xlab = "Mean of normalized counts")

pdf("PDF_figure/CLIP_DESeq_12232019_merged/MAPlot_HFKvHF_peaks.pdf",
    height = 4,
    width = 5)
plotMA(hfk.hf.res,
       main = "MA plot for HFK over HF\nin all peaks (significant over input)",
       xlab = "Mean of normalized counts")
dev.off()
## quartz_off_screen 
##                 2
summary(hfk.hf.res)
## 
## out of 5504 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up)       : 2570, 47%
## LFC < 0 (down)     : 734, 13%
## outliers [1]       : 0, 0%
## low counts [2]     : 0, 0%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
dds_transform <- varianceStabilizingTransformation(dds.peaks.hf.hfk)
rawCountTable_transform <- as.data.frame(assay(dds_transform))
pseudoCount_transform = log2(rawCountTable_transform + 1)
mat.dist = pseudoCount_transform
mat.dist = as.matrix(dist(t(mat.dist)))
mat.dist = mat.dist/max(mat.dist)
setwd("/Users/mizuhi/OneDrive - Harvard University/Haigis Lab/Projects/Halo-Ago2/Halo-Ago-KRas/Raw Data/HEAP-CLIP/HEAP_12232019/Analysis/CLIP_Analysis/Data Visualization")
png('Hierchical_Clustering_merged_peaks.png')
cim(mat.dist, symkey = FALSE, margins = c(6, 6))
suppressMessages(dev.off())
## quartz_off_screen 
##                 2
pdf("PDF_figure/CLIP_DESeq_12232019_merged/Hierchical_Clustering_merged_peaks.pdf")
cim(mat.dist, symkey = FALSE, margins = c(6, 6))
dev.off()
## quartz_off_screen 
##                 2

Final output is following: Hierchical Clustering

plotPCA(dds_transform, intgroup = "tumor", ntop = 500) +   
  geom_text(aes(label=name), vjust = 2, hjust = -0.1) + xlim(-50,60) + ylim(-40,50)

pdf("PDF_figure/CLIP_DESeq_12232019_merged/PCA_merged_peaks.pdf",
    height = 4,
    width = 5)
plotPCA(dds_transform, intgroup = "tumor", ntop = 500) +   
  geom_text(aes(label=name), vjust = 2, hjust = -0.1) + xlim(-50,60) + ylim(-40,50)
dev.off()
## quartz_off_screen 
##                 2
filtered.peaks$hfk.hf.padj <- as.numeric(NA)
filtered.peaks$hfk.hf.log2FC <- as.numeric(NA)
filtered.peaks[rownames(hfk.hf.res), ]$hfk.hf.padj <-
    hfk.hf.res$padj
filtered.peaks[rownames(hfk.hf.res), ]$hfk.hf.log2FC <-
    hfk.hf.res$log2FoldChange

Save the datasets

saveRDS(filtered.peaks, "Datafiles/merged-peaks-filtered-12232019.rds")

Map seed to peaks

mirna.info.family <- readRDS("mirna-info-family-seedmatches.rds")  
assignMirToPeaks <- function(miRNA = mirs, 
                             peaks = es.peaks.utr3,
                             database = mirna.info.family){
  require(BSgenome)
  require(CLIPanalyze)
  require(Biostrings)
  bsgenome <- load.bsgenome("mm10")
  peaks.seq <- get.seqs(bsgenome, peaks)
  peaks$seed.8mer <- as.character(NA)
  peaks$seed.7m8 <- as.character(NA)
  peaks$seed.7A1 <- as.character(NA)
  peaks$seed.6mer <- as.character(NA)
  
  for (i in 1:length(miRNA)){
    mir <- miRNA[i]
    #Prepare seed matches
    mir.6m <- as.character(database[database$miR.family %in% mir, ]$seedmatch.6)
    mir.7m8 <- as.character(database[database$miR.family %in% mir, ]$seedmatch.m8)
    mir.7mA <- as.character(database[database$miR.family %in% mir, ]$seedmatch.A1)
    mir.8m <- as.character(database[database$miR.family %in% mir, ]$seedmatch.8)
    #Filter peaks with 6mer matches
    match <- GRanges(vmatchPattern(mir.6m, peaks.seq))
    #Go back to peaks and map the seed match and extend 1 nt on both directions
    match.extend <- peaks[seqnames(match),]
    match.strand <- as.logical(strand(match.extend) == "+")
    
    start(match.extend[match.strand,]) <- start(match.extend[match.strand,]) + start(match[match.strand,]) - 2
    match.extend[match.strand,] <- resize(match.extend[match.strand,], fix = "start", width = 8)
    
    end(match.extend[!match.strand,]) <- end(match.extend[!match.strand,]) - start(match[!match.strand,]) + 2
    match.extend[!match.strand,] <- resize(match.extend[!match.strand,], fix = "start", width = 8)
    #Assign seed types to these matches
    match.extend.seq <- get.seqs(bsgenome, match.extend)
    match.extend.seq.df <- data.frame(Peaks = names(match.extend.seq),
                                      Sequence = as.character(match.extend.seq))
    
    for (j in 1:nrow(match.extend.seq.df)){
      peak.name <- as.character(match.extend.seq.df[j, "Peaks"])
      seq <- as.character(match.extend.seq.df[j, "Sequence"])
      if (grepl(mir.8m, seq)){
        peaks[peak.name, ]$seed.8mer <- ifelse(is.na(peaks[peak.name, ]$seed.8mer),
                                               mir, 
                                               paste(peaks[peak.name, ]$seed.8mer, mir, sep = ", "))
      } else if (grepl(mir.7m8, seq)){
        peaks[peak.name, ]$seed.7m8 <- ifelse(is.na(peaks[peak.name, ]$seed.7m8),
                                              mir, 
                                              paste(peaks[peak.name, ]$seed.7m8, mir, sep = ", "))
      } else if (grepl(mir.7mA, seq)){
        peaks[peak.name, ]$seed.7A1 <- ifelse(is.na(peaks[peak.name, ]$seed.7A1),
                                              mir, 
                                              paste(peaks[peak.name, ]$seed.7A1, mir, sep = ", "))
      } else {
        peaks[peak.name, ]$seed.6mer <- ifelse(is.na(peaks[peak.name, ]$seed.6mer),
                                                mir, 
                                                paste(peaks[peak.name, ]$seed.6mer, mir, sep = ", "))
      }
    }
  }
  
  return(peaks)
}

miRNAs with mean counts larger than 200 are selected and mapped to peaks

mirna.family.DGE <- readRDS("Datafiles/mirna-counts-deseq-by-family-12232019.rds")
mirs <- subset(mirna.family.DGE, baseMean > 200)
mirs <- mirs$miR.family
#peaks.mirs <- assignMirToPeaks(miRNA = mirs,
#                               peaks = filtered.peaks,
#                               database = mirna.info.family)
#saveRDS(peaks.mirs, "Datafiles/merged-peaks-mirs-200-12232019.rds")
peaks.mirs <- readRDS("Datafiles/merged-peaks-mirs-200-12232019.rds")
targetofmiR <- function(peaks.mir = brain.peaks.mirs,
                        miRNA = "",
                        sitetype = "8mer"){
  peaks.mir.sub <- as.data.frame(peaks.mir[,c("hfk.hf.log2FC", "hfk.hf.padj", 
                                              "seed.8mer", "seed.7m8", "seed.7A1", "seed.6mer")])
  peaks.seedmatch <- lapply(c("seed.8mer", "seed.7m8", "seed.7A1", "seed.6mer"),
                            function(seed){
                              map <- peaks.mir.sub[grepl(miRNA, peaks.mir.sub[,seed]),]
                              map <- rownames(map)
                              map
                            })
  names(peaks.seedmatch) <- c("seed.8mer", "seed.7m8", "seed.7A1", "seed.6mer")

  if (sitetype == "8mer"){
    maps <- peaks.seedmatch[[1]]
  } else if (sitetype == "7mer_above"){
    maps <- unique(unlist(peaks.seedmatch[1:3]))
  } else if (sitetype == "7mer"){
    maps <- unique(unlist(peaks.seedmatch[2:3]))
  } else if (sitetype == "6mer"){
    maps <- peaks.seedmatch[[4]]
  } else {
    print("Please input site type as: 8mer, 7mer_above, 7mer or 6mer")
  }
    return(peaks.mir[maps])
                        }
mirs.peaks <- lapply(mirs,
                     function(mir){
                       targetofmiR(miRNA = mir,
                                   peaks.mir = peaks.mirs,
                                   sitetype = "7mer_above")
                     })
names(mirs.peaks) <- mirs
saveRDS(mirs.peaks, "Datafiles/miRNA-merged-peaks-list-12232019.rds")
lens <- as.data.frame(sapply(mirs.peaks, function(x) length(x)))
mirs.peaks.log2FC.median <- sapply(mirs.peaks,
                            function(list){
                              log2FCs <- list$hfk.hf.log2FC
                              return(median(log2FCs[!is.na(log2FCs)]))
                            })
mirs.peaks.log2FC.mean <- sapply(mirs.peaks,
                            function(list){
                              log2FCs <- list$hfk.hf.log2FC
                              return(mean(log2FCs[!is.na(log2FCs)]))
                            })

mirs.targets.log2FC <- cbind(as.data.frame(mirs.peaks.log2FC.median),
                             as.data.frame(mirs.peaks.log2FC.mean),
                             lens)
colnames(mirs.targets.log2FC) <- c("targets_log2FC_median","targets_log2FC_mean", "N")

mirs.targets.log2FC <- merge(mirs.targets.log2FC, mirna.family.DGE[,c("log2FoldChange", "padj")], by = 0)
colnames(mirs.targets.log2FC)[1] <- "miR.family"

Show any miRNA that has more than 20 peaks matched

df <- subset(mirs.targets.log2FC, N > 20)
p <- ggplot(df, aes(x = log2FoldChange, y = targets_log2FC_median, label = miR.family, size = N)) +
  geom_point(colour = "#1C75BB", alpha = 0.6) +
  #scale_color_manual(fill = "yellow") +
  geom_hline(yintercept = 0, linetype = "dashed", size = 0.3) +
  geom_vline(xintercept = 0, linetype = "dashed", size = 0.3) +
  xlab("miRNA expression Fold change log2 (HFK / HF)") +
  ylab("Median of peak signal changes log2 (HFK / HF)") +
  theme_bw() +
      theme(panel.border = element_blank(),
      panel.background = element_blank(),
      panel.grid.major = element_blank(), 
      panel.grid.minor = element_blank(),
      axis.title.x = element_text(size=12, margin = margin(t = 10)),
      axis.title.y = element_text(size=12, margin = margin(r = 10)),
      axis.text = element_text(size=10),
      axis.line.y = element_line(size = 0.5),
      axis.line.x = element_line(size = 0.5),
      axis.ticks.x = element_line(size = 0),
      axis.ticks.y = element_line(size = 0.5),
      plot.title = element_text(hjust = 0.5, size = 12, face = "bold"))
p1 <- ggplotly(p)
p1
pdf("PDF_figure/CLIP_DESeq_12232019_merged/Median_peak_signal_LFC.pdf",
    width = 5,
    height = 4)
p
dev.off()
## quartz_off_screen 
##                 2
p <- ggplot(df, aes(x = log2FoldChange, y = targets_log2FC_mean, label = miR.family, size = N)) +
  geom_point(colour = "#EC469A", alpha = 0.6) +
  #scale_color_manual(fill = "yellow") +
  geom_hline(yintercept = 0, linetype = "dashed", size = 0.3) +
  geom_vline(xintercept = 0, linetype = "dashed", size = 0.3) +
  xlab("miRNA expression Fold change log2 (HFK / HF)") +
  ylab("Mean of peak signal changes log2 (HFK / HF)") +
  theme_bw() +
      theme(panel.border = element_blank(),
      panel.background = element_blank(),
      panel.grid.major = element_blank(), 
      panel.grid.minor = element_blank(),
      axis.title.x = element_text(size=12, margin = margin(t = 10)),
      axis.title.y = element_text(size=12, margin = margin(r = 10)),
      axis.text = element_text(size=10),
      axis.line.y = element_line(size = 0.5),
      axis.line.x = element_line(size = 0.5),
      axis.ticks.x = element_line(size = 0),
      axis.ticks.y = element_line(size = 0.5),
      plot.title = element_text(hjust = 0.5, size = 12, face = "bold"))
p2 <- ggplotly(p)
p2
pdf("PDF_figure/CLIP_DESeq_12232019_merged/Mean_peak_signal_LFC.pdf",
    width = 5,
    height = 4)
p
dev.off()
## quartz_off_screen 
##                 2

If we only look at the top 40 highly expressed miRNAs

mirna.family.DGE <- mirna.family.DGE[order(-mirna.family.DGE$baseMean),]
mirs.40 <- rownames(mirna.family.DGE)[1:40]

df.40 <- df[df$miR.family %in% mirs.40,]
p <- ggplot(df.40, aes(x = log2FoldChange, y = targets_log2FC_median, label = miR.family, size = N)) +
  geom_point(colour = "#1C75BB", alpha = 0.6) +
  #scale_color_manual(fill = "yellow") +
  geom_hline(yintercept = 0, linetype = "dashed", size = 0.3) +
  geom_vline(xintercept = 0, linetype = "dashed", size = 0.3) +
  xlab("Top 40 miRNA expression Fold change log2 (HFK / HF)") +
  ylab("Median of peak signal changes log2 (HFK / HF)") +
  theme_bw() +
      theme(panel.border = element_blank(),
      panel.background = element_blank(),
      panel.grid.major = element_blank(), 
      panel.grid.minor = element_blank(),
      axis.title.x = element_text(size=12, margin = margin(t = 10)),
      axis.title.y = element_text(size=12, margin = margin(r = 10)),
      axis.text = element_text(size=10),
      axis.line.y = element_line(size = 0.5),
      axis.line.x = element_line(size = 0.5),
      axis.ticks.x = element_line(size = 0),
      axis.ticks.y = element_line(size = 0.5),
      plot.title = element_text(hjust = 0.5, size = 12, face = "bold"))
p3 <- ggplotly(p)
p3
pdf("PDF_figure/CLIP_DESeq_12232019_merged/Median_peak_signal_LFC_top40_miRNA.pdf",
    width = 5,
    height = 4)
p
dev.off()
## quartz_off_screen 
##                 2
p <- ggplot(df.40, aes(x = log2FoldChange, y = targets_log2FC_mean, label = miR.family, size = N)) +
  geom_point(colour = "#EC469A", alpha = 0.6) +
  #scale_color_manual(fill = "yellow") +
  geom_hline(yintercept = 0, linetype = "dashed", size = 0.3) +
  geom_vline(xintercept = 0, linetype = "dashed", size = 0.3) +
  xlab("Top 40 miRNA expression Fold change log2 (HFK / HF)") +
  ylab("Mean of peak signal changes log2 (HFK / HF)") +
  theme_bw() +
      theme(panel.border = element_blank(),
      panel.background = element_blank(),
      panel.grid.major = element_blank(), 
      panel.grid.minor = element_blank(),
      axis.title.x = element_text(size=12, margin = margin(t = 10)),
      axis.title.y = element_text(size=12, margin = margin(r = 10)),
      axis.text = element_text(size=10),
      axis.line.y = element_line(size = 0.5),
      axis.line.x = element_line(size = 0.5),
      axis.ticks.x = element_line(size = 0),
      axis.ticks.y = element_line(size = 0.5),
      plot.title = element_text(hjust = 0.5, size = 12, face = "bold"))
p4 <- ggplotly(p)
p4
pdf("PDF_figure/CLIP_DESeq_12232019_merged/Mean_peak_signal_LFC_top40_miRNA.pdf",
    width = 5,
    height = 4)
p
dev.off()
## quartz_off_screen 
##                 2

Here are the top 10 miRNAs with the most peaks associated.

color.vec <- brewer.pal(name = "Spectral", n = 11)
df <- as.data.frame(filtered.peaks)
df$hfk.hf.logP <- -log10(as.numeric(df$hfk.hf.padj))
len <- sapply(mirs.peaks, function(x) length(x))
mirs.peaks <- mirs.peaks[order(-len)]
mirna <- names(mirs.peaks)
mirs.peaks.names <- lapply(mirs.peaks, function(x) names(x))

for (i in 1:10){
  p <- ggplot() + 
    geom_point(data = df, aes(x = hfk.hf.log2FC, y = hfk.hf.logP), size = 1, alpha = 0.5, color = "grey80") +
    geom_point(data = df[mirs.peaks.names[[mirna[i]]],], aes(x = hfk.hf.log2FC, y = hfk.hf.logP), size = 1, alpha = 0.7, color = "#7570B3") +
    geom_hline(yintercept = -log10(0.05), linetype = "dashed", size = 0.3) +
    geom_vline(xintercept = c(-0.5, 0.5), linetype = "dashed", size = 0.3) +
    xlab("Fold change log2 (HFK / HF)") +
    ylab("-log10(FDR)") +
    theme_bw() +
    theme(panel.border = element_blank(),
      panel.background = element_blank(),
      panel.grid.major = element_blank(), 
      panel.grid.minor = element_blank(),
      axis.title.x = element_text(size=12, margin = margin(t = 10)),
      axis.title.y = element_text(size=12, margin = margin(r = 10)),
      axis.text = element_text(size=10),
      axis.line.y = element_line(size = 0.5),
      axis.line.x = element_line(size = 0.5),
      axis.ticks.x = element_line(size = 0),
      axis.ticks.y = element_line(size = 0.5),
      plot.title = element_text(hjust = 0.5, size = 12, face = "bold")) +
    #guides(fill = guide_legend(title = "miRNA family")) +
    ggtitle(sprintf("%s targets", mirna[i ]))
  print(p)
}
## Warning: Removed 50 rows containing missing values (geom_point).
## Warning: Removed 2 rows containing missing values (geom_point).

## Warning: Removed 50 rows containing missing values (geom_point).
## Warning: Removed 3 rows containing missing values (geom_point).

## Warning: Removed 50 rows containing missing values (geom_point).

## Warning: Removed 50 rows containing missing values (geom_point).

## Warning: Removed 50 rows containing missing values (geom_point).

## Warning: Removed 50 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 50 rows containing missing values (geom_point).

## Warning: Removed 50 rows containing missing values (geom_point).

## Warning: Removed 50 rows containing missing values (geom_point).

## Warning: Removed 50 rows containing missing values (geom_point).

pdf("PDF_figure/CLIP_DESeq_12232019_merged/VolcanoPlot_Peak_by_miRNAFamily.pdf",
    width = 5,
    height = 4)
for (i in 1:10){
  p <- ggplot() + 
    geom_point(data = df, aes(x = hfk.hf.log2FC, y = hfk.hf.logP), size = 1, alpha = 0.5, color = "grey80") +
    geom_point(data = df[mirs.peaks.names[[mirna[i]]],], aes(x = hfk.hf.log2FC, y = hfk.hf.logP), size = 1, alpha = 0.7, color = "#7570B3") +
    geom_hline(yintercept = -log10(0.05), linetype = "dashed", size = 0.3) +
    geom_vline(xintercept = c(-0.5, 0.5), linetype = "dashed", size = 0.3) +
    xlab("Fold change log2 (HFK / HF)") +
    ylab("-log10(FDR)") +
    theme_bw() +
    theme(panel.border = element_blank(),
      panel.background = element_blank(),
      panel.grid.major = element_blank(), 
      panel.grid.minor = element_blank(),
      axis.title.x = element_text(size=12, margin = margin(t = 10)),
      axis.title.y = element_text(size=12, margin = margin(r = 10)),
      axis.text = element_text(size=10),
      axis.line.y = element_line(size = 0.5),
      axis.line.x = element_line(size = 0.5),
      axis.ticks.x = element_line(size = 0),
      axis.ticks.y = element_line(size = 0.5),
      plot.title = element_text(hjust = 0.5, size = 12, face = "bold")) +
    #guides(fill = guide_legend(title = "miRNA family")) +
    ggtitle(sprintf("%s targets", mirna[i ]))
  print(p)
}
## Warning: Removed 50 rows containing missing values (geom_point).
## Warning: Removed 2 rows containing missing values (geom_point).
## Warning: Removed 50 rows containing missing values (geom_point).
## Warning: Removed 3 rows containing missing values (geom_point).
## Warning: Removed 50 rows containing missing values (geom_point).

## Warning: Removed 50 rows containing missing values (geom_point).

## Warning: Removed 50 rows containing missing values (geom_point).

## Warning: Removed 50 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 50 rows containing missing values (geom_point).

## Warning: Removed 50 rows containing missing values (geom_point).

## Warning: Removed 50 rows containing missing values (geom_point).

## Warning: Removed 50 rows containing missing values (geom_point).
dev.off()
## quartz_off_screen 
##                 2

I just want a table of all top 40 miRNAs and the number of peaks that are associated with them with LFC (HFK/HF) > 0 and < 0.

peak.number <- as.data.frame(table(mirs.peaks[[1]]$hfk.hf.log2FC > 0))
peak.number <- t(peak.number[,-1])
colnames(peak.number) <- c("LFC(HFK/HF) < 0", "LFC(HFK/HF) > 0")
for (i in 2: 40) {
  new.peak <- as.data.frame(table(mirs.peaks[[i]]$hfk.hf.log2FC > 0))
  new.peak <- t(new.peak[,-1])
  peak.number <- rbind(peak.number, new.peak)
}

rownames(peak.number) <- names(mirs.peaks)[1:40]
peak.number <- rbind(peak.number, c(sum(df$hfk.hf.log2FC > 0), sum(df$hfk.hf.log2FC < 0)))
rownames(peak.number)[dim(peak.number)[1]] <- "Total"
peak.number
##                                      LFC(HFK/HF) < 0 LFC(HFK/HF) > 0
## let-7-5p/miR-98-5p                               145             382
## miR-200-3p/429-3p                                 67             274
## miR-194-5p                                        64             277
## miR-15-5p/16-5p/195-5p/322-5p/497-5p              83             212
## miR-29-3p                                         57             184
## miR-24-3p                                         56             161
## miR-22-3p                                         61             135
## miR-27-3p                                         43             141
## miR-26-5p                                         47             130
## miR-103-3p/107-3p                                 53             123
## miR-17-5p/20-5p/93-5p/106-5p                      42             132
## miR-196-5p                                        35             119
## miR-141-3p                                        46              97
## miR-23-3p                                         36             100
## miR-484                                           37              79
## miR-30-5p/384-5p                                  52              62
## miR-7-5p                                          17              92
## miR-19-3p                                         33              73
## miR-148-3p/152-3p                                 31              71
## miR-96-5p                                         18              77
## miR-182-5p                                        19              72
## miR-210-3p                                        23              60
## miR-25-3p/32-5p/92-3p/363-3p/367-3p               30              53
## miR-147-3p                                        26              55
## miR-205-5p                                        27              51
## miR-130-3p/301-3p                                 17              61
## miR-194-2-3p/6926-5p/7055-5p                       9              68
## miR-34-5p/449-5p                                  21              54
## miR-185-5p                                        21              52
## miR-378-3p                                        23              50
## miR-30a-3p/30d-3p/30e-3p                          11              58
## miR-671-5p                                        18              51
## miR-192-5p/215-5p                                 24              44
## miR-224-5p                                         7              54
## miR-139-5p                                         7              52
## miR-181-5p                                        13              44
## miR-146ab-5p                                      11              40
## miR-141-5p/6769b-3p                               20              30
## miR-140-5p                                         8              37
## miR-324-3p                                        10              36
## Total                                             NA              NA

SessionInfo

sessionInfo()
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.4
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] plotly_4.9.4                mixOmics_6.12.2            
##  [3] lattice_0.20-44             MASS_7.3-54                
##  [5] RColorBrewer_1.1-2          CLIPanalyze_0.0.10         
##  [7] GenomicAlignments_1.24.0    Rsamtools_2.4.0            
##  [9] GenomicFeatures_1.40.1      AnnotationDbi_1.50.3       
## [11] plyr_1.8.6                  ggrepel_0.9.1              
## [13] gridExtra_2.3               pheatmap_1.0.12            
## [15] Rsubread_2.2.6              Biostrings_2.56.0          
## [17] XVector_0.28.0              DESeq2_1.28.1              
## [19] SummarizedExperiment_1.18.2 DelayedArray_0.14.1        
## [21] matrixStats_0.59.0          Biobase_2.48.0             
## [23] rtracklayer_1.48.0          GenomicRanges_1.40.0       
## [25] GenomeInfoDb_1.24.2         IRanges_2.22.2             
## [27] S4Vectors_0.26.1            BiocGenerics_0.34.0        
## [29] forcats_0.5.1               stringr_1.4.0              
## [31] dplyr_1.0.6                 purrr_0.3.4                
## [33] readr_1.4.0                 tidyr_1.1.3                
## [35] tibble_3.1.2                ggplot2_3.3.3              
## [37] tidyverse_1.3.1             data.table_1.14.0          
## 
## loaded via a namespace (and not attached):
##  [1] colorspace_2.0-1       ellipsis_0.3.2         corpcor_1.6.9         
##  [4] fs_1.5.0               rstudioapi_0.13        farver_2.1.0          
##  [7] bit64_4.0.5            RSpectra_0.16-0        fansi_0.5.0           
## [10] lubridate_1.7.10       xml2_1.3.2             splines_4.0.3         
## [13] cachem_1.0.5           geneplotter_1.66.0     knitr_1.33            
## [16] jsonlite_1.7.2         broom_0.7.7            annotate_1.66.0       
## [19] dbplyr_2.1.1           compiler_4.0.3         httr_1.4.2            
## [22] biosignals_0.1.0       backports_1.2.1        lazyeval_0.2.2        
## [25] assertthat_0.2.1       Matrix_1.3-4           fastmap_1.1.0         
## [28] cli_2.5.0              htmltools_0.5.1.1      prettyunits_1.1.1     
## [31] tools_4.0.3            igraph_1.2.6           gtable_0.3.0          
## [34] glue_1.4.2             GenomeInfoDbData_1.2.3 reshape2_1.4.4        
## [37] rappdirs_0.3.3         Rcpp_1.0.6             cellranger_1.1.0      
## [40] jquerylib_0.1.4        vctrs_0.3.8            crosstalk_1.1.1       
## [43] xfun_0.23              rvest_1.0.0            lifecycle_1.0.0       
## [46] XML_3.99-0.6           zlibbioc_1.34.0        scales_1.1.1          
## [49] hms_1.1.0              curl_4.3.1             yaml_2.2.1            
## [52] memoise_2.0.0          sass_0.4.0             biomaRt_2.44.4        
## [55] stringi_1.6.2          RSQLite_2.2.7          highr_0.9             
## [58] genefilter_1.70.0      BiocParallel_1.22.0    rlang_0.4.11          
## [61] pkgconfig_2.0.3        bitops_1.0-7           evaluate_0.14         
## [64] labeling_0.4.2         htmlwidgets_1.5.3      bit_4.0.4             
## [67] tidyselect_1.1.1       magrittr_2.0.1         R6_2.5.0              
## [70] generics_0.1.0         DBI_1.1.1              pillar_1.6.1          
## [73] haven_2.4.1            withr_2.4.2            survival_3.2-11       
## [76] RCurl_1.98-1.3         modelr_0.1.8           crayon_1.4.1          
## [79] rARPACK_0.11-0         utf8_1.2.1             ellipse_0.4.2         
## [82] BiocFileCache_1.12.1   rmarkdown_2.9          progress_1.2.2        
## [85] locfit_1.5-9.4         grid_4.0.3             readxl_1.3.1          
## [88] blob_1.2.1             reprex_2.0.0           digest_0.6.27         
## [91] xtable_1.8-4           openssl_1.4.4          munsell_0.5.0         
## [94] viridisLite_0.4.0      bslib_0.2.5.1          askpass_1.1